
# Description ------------------------------------------------------------------
#   Goal: Create figures for papers and presentations.
#   Note: Most figures have a "small" version and a normal version. Commented-
#     out size specifications may be relevant to these versions. The "small"
#     version typically fits unrotated on half a vertical sheet of paper, while
#     the non-"small" version (normal?) is meant to fit rotated on a full sheet
#     of paper.

# Setup ------------------------------------------------------------------------
  # Options
  Sys.setenv(TZ = "US/Pacific")
  options(stringsAsFactors = F)
  options(bitmapType = "cairo")
  # Packages
  library(pacman)
  p_load(prism, raster, data.table, lubridate, stringr,
    ggplot2, ggforce, extrafont, Cairo, rmapshaper, ggthemes, viridis,
    magrittr)
  try(ymd(today()))
  # Directories
  dir_project <- "/Users/edwardarubin/Dropbox/Research/MyProjects/NaturalGas/"
  dir_csv     <- paste0(dir_project, "DataCsv/")
  dir_rds     <- paste0(dir_project, "DataR/")
  dir_spatial <- paste0(dir_project, "DataSpatial/")
  # dir_figures <- paste0(dir_project, "Text/ForSubmission/Figures/")
  dir_figures <- paste0(dir_project, "Presentations/Beamer20180131UO/Figures/")
  # Define colors
  very_light_grey <- "#E0E0E0"
  light_grey      <- "#BDBDBD"
  mid_grey        <- "#9E9E9E"
  dark_grey       <- "#616161"
  red_pink        <- "#E40045"
  # purple          <- "#310769"
  purple          <- "#6A1B9A"
  light_purple    <- "#9E17A0"
  # My ggplot themes
  theme_beamer <- theme_bw() + theme(
    panel.border = element_blank(),
    axis.ticks = element_blank(),
    legend.position = 'bottom',
    legend.key = element_blank(),
    panel.grid = element_blank(),
    text = element_text(family = "Roboto", color = "black",
      size = 16),
    title = element_text(size = 18),
    legend.text = element_text(size = 18),
    legend.key.size = unit(2.5, "cm"))
  theme_paper <- theme_bw() + theme(
    panel.border = element_blank(),
    axis.ticks = element_blank(),
    legend.position = 'bottom',
    legend.key = element_blank(),
    panel.grid = element_blank(),
    text = element_text(family = "CharisSIL", color = "black",
      size = 11),
    axis.text = element_text(color = "black"),
    title = element_text(size = 12),
    legend.text = element_text(size = 11),
    legend.key.size = unit(1.5, "cm"))
  # My own save functions
  save_beamer <- function(gg_tmp, name, width = 16 * 0.7, height = 10 * 0.7,
    themes = NULL) {
      ggsave(filename = paste0(name, ".pdf"),
        path = dir_figures,
        plot = gg_tmp + theme_beamer + themes,
        width = width,
        height = height,
        device = cairo_pdf)
      embed_fonts(paste0(dir_figures, name, ".pdf"))
  }
  save_paper <- function(gg_tmp, name, width = 6.4, height = 4,
    themes = NULL) {
      ggsave(filename = paste0(name, ".pdf"),
        path = dir_figures,
        plot = gg_tmp + theme_paper + themes,
        width = width,
        height = height,
        device = cairo_pdf)
      embed_fonts(paste0(dir_figures, name, ".pdf"))
  }

# Figure: Spot, citygate, and baseline prices ----------------------------------
  # Load Henry Hub price data
  henry_dt <- fread(paste0(dir_csv, "naturalGasPricesHenryHub.csv"))
  # Change names
  setnames(henry_dt, c("date", "price"))
  # Format variables
  henry_dt[, `:=`(
    # Convert date to actual date format
    date = mdy(date),
    # Convert price from dollars per million Btu to thm (divide by 10)
    price = price / 10
  )]
  # Add month variable (round date down); drop date
  henry_dt[, month_date := floor_date(date, unit = "month")][, date := NULL]
  # Average price by month
  henry_dt <- henry_dt[, list(price = mean(price, na.rm = T)), month_date]
  # Add type column for type of price
  henry_dt[, type := "henry hub"]
  # Load citygate prices
  citygate_dt <- fread(paste0(dir_csv, "naturalGasPricesCitygate.csv"))
  # Change names
  setnames(citygate_dt, c("date", "price"))
  # Format variables
  citygate_dt[, `:=`(
    # Convert date to actual date format
    date = mdy(date),
    # Convert price from dollars per 1000 cubic feet to thm (divide by 10)
    price = price / 10
  )]
  # Add month variable (round down); drop date
  citygate_dt[, month_date := floor_date(date)][, date := NULL]
  # Add type column for type of price
  citygate_dt[, type := "citygate"]
  # Load the utilities' prices
  utility_dt <- readRDS(paste0(dir_figures, "ggplotPriceSeries.rds"))$data
  # Keep month_date, charge_base, and utility; rename charge_base to price
  utility_dt <- utility_dt[,
    list(month_date, price = charge_base, type = utility)]
  # Stack all price datasets
  price_dt <- rbindlist(list(utility_dt, henry_dt, citygate_dt),
    use.names = T, fill = T)
  # Restrict to 2009-2015 (inclusive)
  price_dt <- price_dt[year(month_date) >= 2009 & year(month_date) <=2015]
  # Plot
  gg_tmp <- ggplot(data = price_dt,
    aes(x = month_date, y = price, color = type, linetype = type)) +
    # geom_path() +
    geom_path(size = 0.4 ) +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = as.numeric(ymd(20090101)),
      color = "grey20", size = 0.4) +
    xlab("Date") + ylab("Price (USD per therm)") +
    ylim(c(0, max(price_dt$price))) +
    scale_color_manual("Price:",
      values = c("#1240AB", "#009999", "grey70", "black"),
      labels = c("Henry Hub spot", "PG&E", "SoCalGas")) +
    scale_linetype_manual("Price:",
      values = c("solid", "dashed", "twodash"),
      labels = c("Henry Hub spot", "PG&E", "SoCalGas"))
  # Save
  save_beamer(gg_tmp, name = "priceComparison")
  # save_paper(gg_tmp, name = "priceComparison", 9, 5.6)
  # save_paper(gg_tmp, name = "priceComparisonSmall",
  #   width = 6.4, height = 3.8,
  #   themes = theme(legend.margin=margin(t=-0.75, r=0.1, b=-0.5, l=0.1, unit="cm")))
  rm(list = c("gg_tmp", "price_dt", "utility_dt", "citygate_dt", "henry_dt"))
  gc()

# Figure: PG&E rate breakdown time series --------------------------------------
  # Load the PG&E rate breakdown data
  pge_dt <- fread(paste0(dir_csv, "pgeRateBreakdown.csv"))
  # Format dates
  pge_dt[, date_eff := mdy(date_eff)]
  # Convert to long form
  pge_dt %<>% melt(id.var = "date_eff", variable.name = "type")
  # Plot the rate breakdown over time
  gg_tmp <- ggplot(data = pge_dt,
    aes(x = date_eff, y = value, color = type, linetype = type)) +
    geom_line() +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = ymd(20000101),
      color = "grey20", size = 0.4) +
    xlab("") +
    ylab("Volumetric charge ($ per therm)") +
    scale_y_continuous(labels = function(x) sprintf("%.2f", x)) +
    scale_color_manual("Price component:",
      values = c(purple, dark_grey, light_grey, "black"),
      labels = c("Procurement", "Baseline transportation",
        "Excess transportation", "PPPS")) +
    scale_linetype_manual("Price component:",
      values = c("solid", "twodash", "twodash", "dashed"),
      labels = c("Procurement", "Baseline transportation",
        "Excess transportation", "PPPS"))
  # Save
  save_beamer(gg_tmp, "pgeTariffSeries",
    themes = theme(
      # legend.direction = "vertical",
      # legend.title = element_text(face = "bold", size = 11),
      legend.title = element_text(face = "bold", size = 16),
      legend.text = element_text(size = 15),
      legend.text.align = 0,
      # axis.text = element_blank(),
      legend.key.size = unit(12, "mm"),
      NULL
    ))
  rm(list = c("gg_tmp", "henry_dt"))
  gc()

# Figure: Henry Hub spot price time series -------------------------------------
  # Load Henry Hub price data
  henry_dt <- fread(paste0(dir_csv, "naturalGasPricesHenryHub.csv"))
  # Change names
  setnames(henry_dt, c("date", "price"))
  # Format variables
  henry_dt[, `:=`(
    # Convert date to actual date format
    date = mdy(date),
    # Convert price from dollars per million Btu to thm (divide by 10)
    price = price / 10
  )]
  # Plot the Henry Hub spot price over time
  gg_tmp <- ggplot(data = henry_dt,
    aes(x = date, y = price)) +
    geom_ribbon(aes(ymin = 0, ymax = price), fill = "grey90") +
    geom_line(color = "grey85", size = 0.3) +
    geom_point(color = "grey25", size = 0.6, shape = 16) +
    # geom_line(color = "grey90", size = 0.1) +
    # geom_point(color = "grey25", size = 0.1, shape = 16) +
    geom_hline(yintercept = 0) +
    # geom_vline(xintercept = min(henry_dt$date) %>% as.numeric()) +
    xlab("Date") +
    ylab("Daily natural gas spot price (USD per therm)") +
    ylim(c(0, max(henry_dt$price)))
  # Save
  save_beamer(gg_tmp, "priceSpot")
  # save_paper(gg_tmp, "priceSpot", 9, 5.6)
  # save_paper(gg_tmp, "priceSpotSmall")
  rm(list = c("gg_tmp", "henry_dt"))
  gc()

# Figure: CA residential natural gas -------------------------------------------
  # Load CA residential natural gas consumption data
  res_dt <- fread(paste0(dir_csv,
    "naturalGasResidentialConsumptionMonthlyCal.csv"))
  setnames(res_dt, c("month", "consumption"))
  # Convert to date
  res_dt[, month := dmy(paste0("1 ", month))]
  gg_tmp <- ggplot(aes(x = month, y = consumption),
    data = res_dt[year(month) >= 2009 & year(month) < 2016]) +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = ymd(20090101),
      color = "grey20", size = 0.4) +
    geom_line(color = dark_grey, alpha = 0.65) +
    geom_point(color = red_pink, alpha = 0.65, size = 1.75, shape = 1) +
    ylim(c(0, max(res_dt[[2]]))) +
    ylab(expression(paste("Monthly natural gas consumption (million ",
      ft^3,")"))) +
    xlab("")
  save_beamer(gg_tmp, "californiaGasConsumption")
  # save_paper(gg_tmp, "californiaGasConsumption", 9, 5.6)
  # save_paper(gg_tmp, "californiaGasConsumptionSmall")
  rm(list = c("res_dt", "gg_tmp"))
  gc()

# Figure: PG&E price regime example --------------------------------------------
  # NOTE: PGE, 2009-01, climate zone 1, had 1.85 therms of daily allowance
  # Load data
  pge_dt <- readRDS(paste0(dir_figures, "ggplotPriceSeries.rds"))$data
  pge_dt <- pge_dt[utility == "pge" & month == 1 & year == 2009]
  # Define price regime structure
  regime_dt <- rbindlist(list(
    data.table(x = 0, xend = 1.85 * 30,
      y = pge_dt$charge_base, yend = pge_dt$charge_base),
    data.table(x = 1.85 * 30, x = 1.85 * 30,
      y = pge_dt$charge_base, yend = pge_dt$charge_excess),
    data.table(x = 1.85 * 30, xend = 200,
      y = pge_dt$charge_excess, yend = pge_dt$charge_excess)
    ))
  # Average price function
  p_avg_fun <- function(x) {
    if(x <= 1.85 * 30) p_avg <- pge_dt$charge_base
    if(x > 1.85 * 30) {
      p_avg <- (1/x) * (pge_dt$charge_base * 1.85 * 30 +
        (x - 1.85 * 30) * pge_dt$charge_excess)
      }
    return(p_avg)
  }
  # The plot
  gg_tmp <- ggplot(data = regime_dt) +
    geom_segment(aes(x = x, xend = xend, y = y, yend = yend,
      color = "1", linetype = "1"),
    size = 0.5) +
    # size = 0.25) +
    geom_line(data = data.table(
      x = seq(0, 200, 0.1),
      y = lapply(X = seq(0, 200, 0.1), FUN = p_avg_fun) %>% unlist()),
      aes(x = x, y = y, color = "2", linetype = "2"),
      size = 0.5) +
      # size = 0.25) +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = 0, color = "grey20", size = 0.4) +
    xlab("Quantity consumed during bill (therms)") +
    ylab("Price (USD per therm)") +
    ylim(c(0, pge_dt$charge_excess * 1.1)) +
    scale_color_manual("",
      values = c("grey70", "black"),
      labels = c("Marginal price", "Average price")) +
    scale_linetype_manual("",
      values = c("solid", "dashed"),
      labels = c("Marginal price", "Average price"))
  # Save
  save_beamer(gg_tmp, "priceRegime")
  # save_paper(gg_tmp, "priceRegime", 9, 5.6)
  # save_paper(gg_tmp, "priceRegimeSmall",
  #   width = 6.4, height = 3.5,
  #   themes = theme(legend.margin=margin(t=-0.5, r=0, b=-0.5, l=0, unit="cm")))
  rm(list = c("pge_dt", "regime_dt", "gg_tmp", "p_avg_fun"))
  gc()

# Figure: PG&E and SoCal price regime series -----------------------------------
  # Load data
  price_dt <- readRDS(paste0(dir_figures, "ggplotPriceSeries.rds"))$data
  # Wide to long
  price_dt <- price_dt[, .(month_date, charge_base, charge_excess, utility)] %>%
    melt(id.var = c("month_date", "utility"),
    variable.name = "price_type", value.name = "price")
  # Create joint variable of utility and price type
  price_dt[, utility_type := paste0(utility, "_", price_type)]
  # The plot
  gg_tmp <- ggplot(data = price_dt,
    aes(x = month_date, y = price, color = utility_type, linetype = utility_type)) +
    geom_line(
      size = 0.5) +
      # size = 0.4) +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = as.numeric(ymd(20090101)),
      color = "grey20", size = 0.4) +
    xlab("Date") +
    ylab("Price (USD per therm)") +
    ylim(c(0, max(price_dt$price))) +
    scale_color_manual("Price:",
      values = rep(c("#6A1B9A", "#FC845F"), each = 2),
      labels = c("PG&E: Baseline", "PG&E: Excess", "SoCalGas: Baseline", "SoCalGas: Excess")) +
    scale_linetype_manual("Price:",
      values = rep(c("solid", "twodash"), 2),
      labels = c("PG&E: Baseline", "PG&E: Excess", "SoCalGas: Baseline", "SoCalGas: Excess")) +
    theme(legend.box = "horizontal")
  # Save
  save_beamer(gg_tmp, "priceRegimeSeries",
  # save_paper(gg_tmp, "priceRegimeSeries", 9, 5.6)
  # save_paper(gg_tmp, "priceRegimeSeriesSmall",
  #   width = 6.4, height = 3.8,
  #   themes = theme(legend.margin=margin(t=-0.75, r=0.1, b=-0.5, l=0.1, unit="cm")))
    themes = theme(
      legend.title = element_text(face = "bold"),
      legend.text = element_text(size = 16),
      legend.key.width = unit(1.5, "cm"),
      legend.key.height = unit(2.5, "cm"),
      NULL))
  rm(list = c("price_dt", "gg_tmp"))
  gc()

# Figure: Allowance series -----------------------------------------------------
  # NOTE: PG&E "R" zone and SoCalGas "1" zone
  # Load data
  allow_dt <- readRDS(paste0(dir_figures, "ggplotAllowanceSeries.rds"))$data
  # Add end of month
  allow_dt[, month_end :=
    ceiling_date(month_date, unit = "month", change_on_boundary = T) - 1]
  # Copy the month's allowance (horizonal line within month)
  allow_dt[, allow_end := allow_ind]
  # Copy dataset for movement between months
  copy_dt <- copy(allow_dt)
  copy_dt[, month_date := month_end]
  setorder(copy_dt, utility, month_date)
  copy_dt[, allow_end := shift(allow_ind, 1L, NA, "lead"), by = utility]
  # Bind datasets
  allow_dt <- rbindlist(list(allow_dt, copy_dt)); rm(copy_dt); gc()
  setorder(allow_dt, utility, month_date)
  # Plot
  gg_tmp <- ggplot(data = allow_dt,
    aes(color = utility, linetype = utility)) +
    geom_segment(aes(x = month_date, xend = month_end,
      y = allow_ind, yend = allow_end),
      # size = 0.5) +
      size = 0.4) +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = as.numeric(ymd(20090101)),
      color = "grey20", size = 0.4) +
    xlab("") +
    ylab("Daily allowance (therms)") +
    ylim(c(0, max(allow_dt$allow_ind))) +
    scale_color_manual("",
      values = c("#6A1B9A", "#FC845F"),
      labels = c("PG&E", "SoCalGas")) +
    scale_linetype_manual("",
      values = c("longdash", "twodash"),
      labels = c("PG&E", "SoCalGas"))
  # Save
  save_beamer(gg_tmp, "allowanceSeries")
  # save_paper(gg_tmp, "allowanceSeries", 9, 5.6)
  embed_fonts(paste0(dir_figures, "allowanceSeries.pdf"))
  # save_paper(gg_tmp, "allowanceSeriesSmall",
  #   width = 6.4, height = 3.8,
  #   themes = theme(legend.margin=margin(t=-1, r=0, b=-0.5, l=0, unit="cm")))
  # embed_fonts(paste0(dir_figures, "allowanceSeriesSmall.pdf"))
  rm(gg_tmp); gc()

# Figure: Climate zone map -----------------------------------------------------
  # Load additional packages
  library(rgdal)
  library(broom)
  # Load shapefile of climate zones
  ca_shp <- readOGR(
    dsn = paste0(dir_spatial, "CA_Building_Standards_Climate_Zones"),
    layer = "CA_Building_Standards_Climate_Zones")
  ca_shp %<>% ms_simplify(keep = 0.1, keep_shapes = T)
  # Tidy to plot
  ca_dt <- tidy(ca_shp, group = "Zone")
  # Plot
  gg_tmp <- ggplot(data = ca_dt, aes(x = long, y = lat, group = group)) +
    geom_path(size = 0.2) +
    # geom_path(size = 0.1) +
    xlab("") + ylab("") +
    coord_quickmap()
  # Save
  save_beamer(gg_tmp, "climateZones",
    themes = theme(axis.text = element_blank()))
  # save_paper(gg_tmp, "climateZones", 6, 8,
  # save_paper(gg_tmp, "climateZonesSmall", 3, 4,
  #   themes = theme(axis.text = element_blank()))
  rm(list = c("ca_shp", "gg_tmp"))

# Figure: Common zip code map --------------------------------------------------
  # Load additional packages
  p_load(rgdal, broom, rgeos, maptools)
  # Load the zip code shapefile
  zip_shp <- readOGR(
    dsn = paste0(dir_spatial, "CaliforniaZipCodes"),
    layer = "CaliforniaZipCodes")
  # Load the common cities CSV
  common_dt <- fread(paste0(dir_csv, "commonCityZips.csv"))
  # Convert zips to character
  common_dt[, zip := as.character(zip)]
  # Take subset of shapefile for zip codes that are in the common cities file
  zip_shp %<>% subset(ZIP_CODE %in% common_dt$zip)
  # Convert to data.table
  zip_dt <- tidy(zip_shp, region = "ZIP_CODE") %>% data.table()
  # Join information
  zip_dt %<>% merge(y = common_dt, by.x = "id", by.y = "zip",
    all.x = T, all.y = F, sort = F)
  # Add column combining utilities' presence
  zip_dt[pge == F & socal == F, utility := "neither"]
  zip_dt[pge == T & socal == F, utility := "pge"]
  zip_dt[pge == F & socal == T, utility := "socal"]
  zip_dt[pge == T & socal == T, utility := "both"]
  # Cities' centers
  city_centers <- zip_dt[, list(long = mean(long), lat = mean(lat)), by = city]
  # More adjustements
  city_centers[city == "Del Rey", lat := lat + 0.03]
  city_centers[city == "Fellows", lat := lat + 0.10]
  city_centers[city == "Tehachapi", lat := lat - 0.03]
  city_centers[city == "Paso Robles", lat := lat + 0.01]
  city_centers[city == "Templeton", lat := lat - 0.1]
  # Drop "neither" and order
  zip_dt <- zip_dt[utility != "neither"]
  zip_dt[, utility_fac := factor(utility,
    levels = c("pge", "both"),
    ordered = T,
    labels = c("Only PG&E", "Both PG&E and SoCalGas"))]
  # Plot
  gg_tmp <- ggplot(data = zip_dt,
    aes(x = long, y = lat, group = group)) +
    geom_polygon(aes(fill = utility_fac), alpha = 0.9) +
    geom_path(color = "grey93", size = 0.15, alpha = 0.9) +
    xlab("") + ylab("") +
    scale_fill_manual("Utility presence:",
    # scale_fill_manual("Utilities serving the zip code:",
      # values = c("#4DB6AC", "#B0BEC5"),
      # values = c("#B0BEC5", "#4DB6AC"),
      values = c("#ffc107", "#e40045")) +
    geom_text(data = city_centers, family = "Roboto", size = 3.6,
    # geom_text(data = city_centers, family = "CharisSIL", size = 3.5,
      aes(x = long, y = lat, label = city, group = city)) +
    xlim(zip_dt$long %>% range()) +
    ylim(zip_dt$lat %>% range()) +
    coord_quickmap()
  # Save
  # save_beamer(gg_tmp, "studyAreaNoLegend", width = 14 * 0.7, height = 14 * 0.7,
  save_beamer(gg_tmp, "studyAreaWithoutLegend", width = 6.1, height = 6.1,
    themes = theme(
      legend.position = "none",
      axis.text = element_blank(),
      legend.key.size = unit(1, "cm"),
      # legend.title = element_text(face = "bold", size = 12),
      # legend.text = element_text(size = 11),
      NULL
    ))
  # save_paper(gg_tmp, "studyArea", width = 6, height = 6,
  #   themes = theme(
  #     legend.position = "none",
  #     axis.text = element_blank(),
  #     legend.key.size = unit(1, "cm"),
  #     # # Add panel to match highlighted rectangle from service area map (below)
  #     # panel.border  = element_rect(
  #     #   fill = NA, color = "grey70", size = 0.25, linetype = "dashed"),
  #     NULL
  #   ))
  # rm(list = c("zip_shp", "zip_dt", "common_dt", "city_centers", "gg_tmp"))
  # gc()

# Figure: Extended zip code map ------------------------------------------------
  # Load additional packages
  p_load(rgdal, broom, rgeos, maptools)
  # Load the neighbors dataset
  nbr_dt <- readRDS(paste0(dir_rds, "zipCodeGroups.rds"))
  # Keep the desired groups (0-3 currently)
  nbr_dt <- nbr_dt[group < 4]
  setorder(nbr_dt, group, zip)
  nbr_dt <- nbr_dt[, .(zip, nbr_grp = group)]
  # Load the zip code shapefile
  zip_shp <- readOGR(
    dsn = paste0(dir_spatial, "CaliforniaZipCodes"),
    layer = "CaliforniaZipCodes")
  # Subset to the zip codes in our groups
  zip_shp <- subset(zip_shp, ZIP_CODE %in% nbr_dt[,zip])
  # Convert to data.table
  zip_dt <- tidy(zip_shp, region = "ZIP_CODE") %>% data.table()
  # Join information
  zip_dt %<>% merge(y = nbr_dt, by.x = "id", by.y = "zip",
    all.x = T, all.y = F, sort = F)
  # Set order
  zip_dt[, nbr_fac := factor(nbr_grp,
    levels = 0:3,
    ordered = T,
    labels = c(
      "Common Zips: zip codes served by both utilities",
      "Neighbors 1: 'Common Zips' and their neighbors",
      "Neighbors 2: 'Neighbors 1' and their neighbors",
      "Neighbors 3: 'Neighbors 2' and their neighbors"))]
  # Plot
  gg_tmp <- ggplot(data = zip_dt,
    aes(x = long, y = lat, group = group)) +
    geom_polygon(aes(fill = nbr_fac), alpha = 0.9) +
    # geom_path(color = "grey90", size = 0.1, alpha = 0.9) +
    geom_path(color = "white", size = 0.1, alpha = 0.9) +
    xlab("") + ylab("") +
    scale_fill_manual("Zip-code group",
      values = inferno(4, direction = 1, end = 0.9),
      labels = expression(
        paste(bold("Common Zips: "), "zip codes served by both utilities"),
        paste(bold("Neighbors 1: "), italic("Common Zips "), "and their neighbors"),
        paste(bold("Neighbors 2: "), italic("Neighbors 1 "), "and their neighbors"),
        paste(bold("Neighbors 3: "), italic("Neighbors 2 "), "and their neighbors")
      )) +
    xlim(zip_dt$long %>% range()) +
    ylim(zip_dt$lat %>% range()) +
    coord_quickmap()
  # Save
  save_beamer(gg_tmp, "neighbors", width = 14 * 0.7, height = 14 * 0.7,
  # save_paper(gg_tmp, "neighbors", width = 6.5, height = 6.5,
    themes = theme(
      legend.direction = "vertical",
      # legend.title = element_text(face = "bold", size = 11),
      legend.title = element_text(size = 18),
      legend.text.align = 0,
      axis.text = element_blank(),
      legend.key.size = unit(9, "mm"),
      NULL
    ))
  rm(list = c("zip_shp"))
  gc()

# Figure: Service areas map ----------------------------------------------------
  # Load additional packages
  p_load(rgdal, broom, rgeos, maptools, viridis, rmapshaper)
  # Load the zip code shapefile
  zip_shp <- readOGR(
    dsn = paste0(dir_spatial, "CaliforniaZipCodes"),
    layer = "CaliforniaZipCodes")
  # Take union of the zip codes
  ca_shp <- gUnaryUnion(zip_shp)
  # Load the zip code and utility list
  utility_dt <- fread(paste0(dir_csv, "utilityZipCodeMap.csv"))
  # Take subsets for each set of utilities and then take union
  pge_shp <- subset(zip_shp,
    ZIP_CODE %in% utility_dt[pge == T & socal == F, zip]) %>% gUnaryUnion()
  socal_shp <- subset(zip_shp,
    ZIP_CODE %in% utility_dt[pge == F & socal == T, zip]) %>% gUnaryUnion()
  both_shp <- subset(zip_shp,
    ZIP_CODE %in% utility_dt[pge == T & socal == T, zip]) %>% gUnaryUnion()
  # Reduce the resolution of each shapefile
  pge_simple <- ms_simplify(pge_shp, keep = 0.01)
  socal_simple <- ms_simplify(socal_shp, keep = 0.01)
  both_simple <- ms_simplify(both_shp, keep = 0.01)
  ca_simple <- ms_simplify(ca_shp, keep = 0.01, keep_shapes = T)
  # Convert shapefiles to data tables
  pge_dt <- pge_shp %>% tidy(region = "ZIP_CODE") %>% data.table()
  socal_dt <- socal_shp %>% tidy(region = "ZIP_CODE") %>% data.table()
  both_dt <- both_shp %>% tidy(region = "ZIP_CODE") %>% data.table()
  ca_dt <- ca_shp %>% tidy(region = "ZIP_CODE") %>% data.table()
  pge_sdt <- pge_simple %>% tidy(region = "ZIP_CODE") %>% data.table()
  socal_sdt <- socal_simple %>% tidy(region = "ZIP_CODE") %>% data.table()
  both_sdt <- both_simple %>% tidy(region = "ZIP_CODE") %>% data.table()
  ca_sdt <- ca_simple %>% tidy(region = "ZIP_CODE") %>% data.table()
  # Plot
  gg_tmp <- ggplot(data = pge_sdt, aes(x = long, y = lat, group = group)) +
    geom_polygon(data = pge_sdt, aes(fill = "1")) +
    geom_polygon(data = socal_sdt, aes(fill = "2")) +
    geom_polygon(data = both_sdt, aes(fill = "3")) +
    geom_path(data = ca_sdt, size = 0.1, color = "grey20") +
    # # Box highlighting the study area
    # annotate("rect",
    #   xmin = both_shp@bbox["x","min"],
    #   xmax = both_shp@bbox["x","max"],
    #   ymin = both_shp@bbox["y","min"],
    #   ymax = both_shp@bbox["y","max"],
    #   fill = NA, color = "grey70", size = 0.25, linetype = "dashed") +
    xlab("") + ylab("") +
    scale_fill_manual("Utility presence:",
      values = c("#ffc107", "#310769", "#e40045"),
      labels = c("PG&E", "SoCalGas", "PG&E and SoCalGas")) +
    coord_quickmap()
  # Save
  save_beamer(gg_tmp, "serviceAreas", width = 14 * 0.7, height = 14 * 0.7,
  # save_paper(gg_tmp, "serviceAreas", width = 6, height = 7,
    themes = theme(
      axis.text = element_blank(),
      legend.key.size = unit(1, "cm"),
      legend.title = element_text(face = "bold"),
      NULL
    ))
  # ggsave(filename = "serviceAreas.png", plot = gg_tmp, path = dir_figures,
  #   width = 14 * 0.7, height = 14 * 0.7)
  rm(list = c("zip_shp", "utility_dt", "zip_dt", "common_dt", "gg_tmp"))
  gc()

# Calendar example -------------------------------------------------------------
  # Create dataset for calendar
  # Days in weeks
  cal_dt <- expand.grid(x = 1:7, y = 1:20) %>% data.table()
  # Days of months
  cal_dt[, date := c(1:31, 1:30)]
  # Months
  cal_dt[, month :=
    mapply(
      FUN = function(i,j) rep(i, each = j),
      i = 0:9,
      j = rep(c(31,30), 5),
      SIMPLIFY = F) %>% unlist()]
  # Billing periods (start on the 6th of the month)
  cal_dt[, bill := c(rep(0, each = 10),
    mapply(
      FUN = function(i,j) rep(i, each = j),
      i = 1:10,
      j = rep(c(31,30), 5),
      SIMPLIFY = F) %>% unlist())
      ]
  # # Keep the first 4 "bills"
  # cal_dt <- cal_dt[bill < 5]
  # # HACK: Fill out the rest of the last week (currently only needs one more day)
  # cal_dt <- rbindlist(list(cal_dt,
  #   data.table(x = 7, y = 19, date = 11, month = 4, bill = 5)
  # ))
  cal_dt[, bill_f := factor(bill,
    levels = c(0:5),
    labels = c("Lag 4", "Lag 3", "Lag 2", "Lag 1", "Current", "Lead 1"),
    ordered = T)]
  # Update 'y' and create a 'y' that takes into account month
  cal_dt[, y := -y]
  cal_dt[, y_m := y - month]
  # Plot
  gg_tmp <- ggplot(data = cal_dt, aes(x,y_m)) +
    # Main grid
    geom_tile(width = 0.94, height = 0.96,
      aes(fill = bill_f), color = "white") +
    # Labels
    geom_tile(width = 0.96, height = 0.96, data = cal_dt[bill == 4 & date == 17,],
      aes(x,y_m), fill = NA, color = "black", size = 0.9) +
    geom_tile(width = 0.96, height = 0.96, data = cal_dt[bill == 4 & date == 1,],
      aes(x,y_m), fill = NA, color = "black", size = 0.9) +
    # Calendar text
    geom_text(aes(label = date), color = c(
      rep("grey70", 10),
      rep("white", 31),
      rep("white", 30),
      rep("white", 31),
      rep("black", 30),
      rep("grey70", 8),
      NULL),
     # family = "Charter", size = 3.25) +
     family = "Roboto", size = 3.15) +
    # Label "Bill received"
    geom_segment(x = 4.5, xend = 8, y = -18.508, yend = -18.508,
      color = "black", size = 0.3) +
    annotate(geom = "Text", x = 8.1, y = -18.5, color = "black", hjust = 0,
      # label = "Bill received", family = "Charter") +
      label = "Bill received", family = "Roboto") +
    # Label "Bill due"
    geom_segment(x = 4.5, xend = 8, y = -21.508, yend = -21.508,
      color = "black", size = 0.3) +
    annotate(geom = "Text", x = 8.1, y = -21.5, color = "black", hjust = 0,
      # label = "Bill due",  family = "Charter") +
      label = "Bill due",  family = "Roboto") +
    # # Highlight days with circles
    # geom_point(data = cal_dt[bill == 4 & date == 17,],
    #   color = "#6A1B9A", size = 12) +
    # geom_text(data = cal_dt[bill == 4 & date == 17,], aes(label = date),
    #   color = "white", size = 4) +
    # geom_point(data = cal_dt[bill == 4 & date == 1,],
    #   color = "#FF5722", size = 12) +
    # geom_text(data = cal_dt[bill == 4 & date == 1,], aes(label = date),
    #   color = "white", size = 4) +
    # Days of week
    annotate(geom = "text", x = 1:7, y = 0, color = "black",
      # label = c("S", "M", "T", "W", "T", "F", "S"), family = "CharisSIL") +
      label = c("S", "M", "T", "W", "T", "F", "S"), family = "Roboto") +
    # Line below the days of week
    annotate(geom = "segment", x = 0.55, xend = 7.45, y = -0.5, yend = -0.5,
      color = "grey75", size = 0.5) +
    # Limits
    xlim(0,11) +
    # Aesthetics
    coord_equal() +
    labs(x = "", y = "") +
    scale_fill_manual(
      "Billing period",
      breaks = c("Lag 3", "Lag 2", "Lag 1", "Current"),
      values = c("white", "grey15", "grey40", "grey65", "grey85", "white")) +
    theme_bw() + theme(
      panel.border = element_blank(),
      axis.ticks = element_blank(),
      legend.key = element_rect(color = "white", size = 1.2),
      panel.grid = element_blank(),
      # text = element_text(family = "CharisSIL", color = "black", size = 11),
      text = element_text(family = "Roboto", color = "black", size = 11),
      axis.text = element_blank(),
      title = element_text(size = 11),
      # legend.text = element_text(size = 11, family = "Charter", face = "plain"),
      # legend.title = element_text(size = 11, family = "Charter", face = "italic"),
      legend.text = element_text(size = 11, family = "Roboto", face = "plain"),
      legend.title = element_text(size = 11, family = "Roboto", face = "bold"),
      legend.key.size = unit(0.5, "cm"))
    # Save for paper
    # ggsave(
    #   plot = gg_tmp,
    #   path = dir_figures,
    #   filename = "calendar.pdf",
    #   width = 6,
    #   height = 7.25,
    #   device = cairo_pdf)
    # Save for beamer
    ggsave(
      plot = gg_tmp,
      path = dir_figures,
      filename = "calendar.pdf",
      width = 6,
      height = 7.25,
      device = cairo_pdf)

# NOTE: Used QGIS for this map. ggplot2 could not project with raster and
# geom_tile() never finished running.
# # Figure: PRISM daily data map -------------------------------------------------
#   # Choose folder for PRISM
#   options(prism.path = paste0(dir_spatial, "/PRISM"))
#   # Download a (tax) day of temperature from PRISM
#   get_prism_dailys(
#     type = "tmean",
#     minDate = "2010-06-15",
#     maxDate = "2010-06-15",
#     keepZip = F)
#   # Load the PRISM data
#   temp_r <- raster(str_subset(ls_prism_data(absPath = T)[,2], "20100615"))
#   # Convert the raster to a data table
#   temp_dt <- rasterToPoints(temp_r) %>% data.table()
#   # Rename
#   setnames(temp_dt, c("x", "y", "temp_c"))
#   # Convert temperature from C to F
#   temp_dt[, temp_f := temp_c * 1.8  + 32]
#   # Map/plot it
#   gg_tmp <- ggplot(data = temp_dt, aes(x, y, fill = temp_f)) +
#     geom_raster() +
#     scale_fill_viridis("Mean temperature (F)", option = "B") +
#     coord_map(projection = "lambert", lat0 = 20, lat1 = 60) +
#     theme_map()
#   save_paper(gg_tmp, "mapPRISM", width = 8, height = 6,
#     themes = theme(
#       legend.position = "none",
#       axis.text = element_blank(),
#       legend.key.size = unit(1, "cm"),
#       NULL
#     ))
